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The renormalization of the band structure at zero temperature due to electron-phonon coupling 
is explored in diamond, BN, LiF and MgO crystals. We implement a dynamical scheme to com¬ 
pute the frequency-dependent self-energy and the resulting quasiparticle electronic structure. Our 
calculations reveal the presence of a satellite band below the Fermi level of LiF and MgO. We show 
that the renormalization factor (Z), which is neglected in the adiabatic approximation, can reduce 
the zero-point renormalization (ZPR) by as much as 40%. Anharmonic effects in the renormalized 
eigenvalues at finite atomic displacements are explored with the frozen-phonon method. We use a 
non-perturbative expression for the ZPR, going beyond the Allen-Heine-Cardona theory. Our results 
indicate that high-order electron-phonon coupling terms contribute significantly to the zero-point 
renormalization for certain materials. 
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The electron-phonon coupling is at the heart of nu¬ 
merous phenomena such as optical absorption^’^, ther¬ 
moelectric transport^, and superconductivityIt is 
also a crucial ingredient in basic electronic structure 
calculations, giving renormalized quasiparticle energies 
and lifetimes. This renormalization causes the temper¬ 
ature dependence of the band gap of semiconductors®, 
and accounts for the zero-point renormalization (ZPR), 
while the lifetime broadenings are observed through the 
electron mobility®’and in photo-absorption/emission 
experiments^^. 

Obtaining the quasiparticle structure from first prin¬ 
ciples has been a challenge, addressed for the first 
time for bulk silicon by King-Smith et al}'^ , in 1989, 
using density functional theory (DFT), with a mixed 
frozen-phonon supercell and linear response approach. 
These authors pointed the inadequate convergence of 
their results with respect to phonon wavevector sam¬ 
pling, due to the limited available computing capabil¬ 
ities. Fifteen year passed, before Capaz et al. com¬ 
puted it for carbon nanotubes^® using DFT with frozen 
phonons. At variance with the frozen-phonon approach, 
the theory of Allen, Heine and Gardona (AHG)^^“^® 
casts the renormalization and the broadening in terms 
of the first-order derivatives of the effective poten¬ 
tial with respect to atomic positions. Used initially 
with empirical potentials, tight-binding or semi-empirical 
pseudopotentials^^^®®, AHG was then implemented with 
the density functional perturbation theory (DFPT)®^“®^, 
providing an efficient way to compute the phonon band 
structure and the electron-phonon coupling altogether. 
This powerful technique allowed A. Marini to com¬ 
pute, from first principles, temperature-dependent opti¬ 
cal properties®®. 

While DFPT has been widely applied to predict struc¬ 
tural and thermodynamical properties of solids®®, few 
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studies have used it to compute the phonon-induced 
renormalization of the band structure. The scarcity of 
experimental data is at least partly responsible for this 
imbalance. Whereas the phonon spectrum is commonly 
measured through Raman spectroscopy and neutron¬ 
scattering experiments, evaluating the ZPR requires low- 
temperature ellipsometry measurements or isotope sub¬ 
stitutions, which are less abundant in literature. From 
a theoretical point of view, the calculation of the ZPR 
relies on several assumptions that we will be addressing 
in this paper. 

We identify two kinds of approximations. The first 
kind are those regarding the treatment of the electron- 
electron interactions, which is achieved in DFT through 
the Hartree and the exchange-correlation potentials. It 
was shown that the strength of the electron-phonon inter¬ 
action was highly sensitive to the choice of the exchange- 
correlation functional®^. Subsequent GW calculations 
confirmed that standard functionals such as the local- 
density approximation (LDA) tend to underestimate the 
electron-phonon coupling by as much as 30%®®“®^. 

The second kind of approximations are those made on 
the self-energy of the electron-phonon interaction. One, 
for example, usually performs the rigid-ion approxima¬ 
tion, assuming that the second-order derivative of the 
Hamiltonian is diagonal in atom sites. This approxima¬ 
tion proved to be valid in the case of crystals®^’®®, but 
notably fails for diatomic molecules®®. 

Another assumption is the adiabatic approximation, 
which implies that the phonon population can be treated 
as a static perturbation. One would typically compute 
the real part of the self-energy in a static way, and use 
a dynamical expression to compute the imaginary part 
and obtain the electronic lifetimes®^. The adiabatic ap¬ 
proximation breaks downs in certain materials such as 
diamond and polyacetylene, as pointed out by Gannuc- 
cia and Marini®®’®®. By considering the frequency depen¬ 
dence of the self-energy, they showed that the electron- 
phonon interaction smears out the energy levels, even 
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obliterating the band structure. 

Finally, the harmonic approximation is the assump¬ 
tion that the total energy and electronic eigenvalues 
vary quadratically with atomic displacements, which jus¬ 
tifies the use of a second-order perturbation theory. 
Higher order expansions have been used to compute 
phonon wavefunctions, energies, and thermal expansion 
coefficients^^^^®, but its impact on the ZPR was rarely 
investigated. 

In this work, we compute the ZPR and the quasipar¬ 
ticle lifetimes of the band structure of diamond, BN, 
LiF and MgO. We show that the inclusion of dynam¬ 
ical effects in the AHC theory is important to obtain 
correct quasiparticle energies and broadenings. We also 
study the impact of anharmonic effects in the electronic 
energies by means of frozen-phonon calculations, and 
show that high-order terms do contribute to the electron- 
phonon coupling in certain cases. 

All calculations are performed with the Abinit code"*®. 
For simplicity, we use an LDA exchange-correlation func¬ 
tional. We do not expect this approach to fully capture 
the strength of the coupling, as does the GW method. 
Rather, it allows us to evaluate the impact of several com¬ 
monly adopted approximations to the electron-phonon 
coupling self-energy. 


I. DYNAMICAL DFPT 


The dynamical AHC theory is derived by expanding 
the starting Hamiltonian Hq up to second order in atomic 
displacements as (using atomic units) 
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where c\ and c\ are electron creation and destruction 
operators, and A^ = + a^_^, such that A^j rep¬ 

resents a phonon displacement operator. The electronic 
states A with wavefunctions '0 a and energies e° are spec¬ 
ified by a wavevector Ica, a band index 6 a, and spin cta, 
while the phonon modes v with frequencies are spec¬ 
ified by a wavevector and a branch index mj^, and A 
is the number of phonon wavevectors. 

The first-order perturbation potential is formed with 
derivatives of the Hamiltonian with respect to atomic 
displacements along a particular phonon mode as 

Rji) =V.Ho = J2 Ho, (2) 

iKj 


where I labels a unit cell with lattice vector R/, k an atom 
within a unit cell, j a cartesian direction, and is the 
phonon displacement vector. The second-order pertur¬ 
bation potential is then V^} (r) = iVj,V*,Ao. 

Within the DFT and DFPT approaches of this work, 
Hq is an electron-only Hamiltonian, and the phonon 
perturbations involve the derivatives of the local self- 
consistent potential. Since the electronic density re¬ 
sponds to the atomic displacements by screening the ions’ 
potential, these perturbations must be evaluated self- 
consistently. Furthermore, the phonon displacement vec¬ 
tors are a priori unknown and must be computed by di- 
agonalization of the the dynamical matrix. Other formu¬ 
lations of Ho could include the many-body interactions 
between the electrons and the ions^^ . These alternatives 
however fall beyond the scope of this work. 

Following the usual many-body perturbation theory 
the electron-phonon self-energy at second order is the 
sum of the Fan and the Debye-Waller terms: 

Sf,,(a;) = . (3) 

The dynamical Fan term is given by 


= E ^ E 1 ^^") 

V A" 


where and /a are boson and fermion occupation fac¬ 
tors, and 77 is a small parameter which is real and pos¬ 
itive. This parameter maintains causality by giving the 
correct sign to the imaginary part of the quasiparticle 
energies. It also smooths the frequency dependence of 
the self-energy when a finite sampling of phonon modes 
is used (see appendix A). We note that the periodicity of 
the phonon perturbation potential restricts the summa¬ 
tion over intermediate states to those at the k-point given 


nyT) + fy,{T) ^ n,{T) + l- fy,{T) 1 

— £^„ + LOi, + ir] sgn(w) uj — e^„ — + ip sgn(a;)J ’ 

( 4 ) 


by kA" -I- = kA = kA'. In Eq. (4) and in the remaining 

of this work, all the summations over the phonon modes 
are implicitly normalized by the number of wavevectors 
used to sample the Brillouin zone. 

The frequency-independent Debye-Waller term is for¬ 
mally defined as 

= E ^ [2n.(r) + 1 ], (5) 
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which also implies = k^'. Within the rigid-ion ap¬ 
proximation, the Debye-Waller term can be computed 
using only the matrix elements of in a form similar 
to the Fan term, thanks to translational invariance^^. 

The interacting Green’s function is the solution of the 
Dyson equation involving the full electron-phonon self¬ 
energy, which is diagonal in (k^jk^')- If the bands are 
well separated in energy, then the Green’s function can 
be approximated with the diagonal elements of the self¬ 
energy as 


GA(a;)« (a;-e^-Ef(a;)) \ (6) 

where we use the shorthand By con¬ 

sidering only the diagonal elements of the self-energy, we 
disregard the possible change in the electronic density 
resulting from the electron-phonon interaction. Taking 
this effect into account would involve solving the Dyson 
equation for the Green’s function to obtain a new den¬ 
sity, and applying the change in the DFT self-consistent 
potential perturbatively. We are assuming however that 
the change in the one-electron state densities among the 
set of occupied bands (and among the set of unoccupied 
bands) would compensate, and that the additional per¬ 
turbative terms (that is, excitations accross the gap) be 
negligible. Besides, we stress that the Green’s function 
only needs to be corrected once, since this procedure is 
a second-order perturbative approach. Any attempt for 
self-consistency in the calculation of the Green’s function 
and the self-energy belongs to a higher-order treatment. 

From the imaginary part of the Green’s function, one 
obtains the spectral function 

. . ^ 1 PmSf(a;)| 

7r[a;-£0-ineEf(w)]2 + amEf(cu)2’ ^ 

which directly relates to the signal observed in ARPES 
experiments. The quasiparticle energies £\ are defined as 
the position of the principal peak of ^>( 0 ;). Neglecting 
the frequency dependence of the imaginary part of the 
self-energy, the maximum of the spectral function is at 
the energy given by 

£a = 4+5IeEf (£a)- (8) 


Assuming furthermore that the quasiparticle energies are 
close to the bare electronic energies, the latter can be 
corrected perturbatively as 




where 
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du 
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is the renormalization factor. This procedure accounts 
for a linearization of the self-energy near the bare eigen¬ 
value, as illustrated in Fig 1. 



FIG. 1. Upper: Real and imaginary part of the self-energy 
for the top of the valence bands (VB) of LiF. The vertical 
line indicates the bare eigenvalue, the x = y line gives the 
renormalized eigenvalue at the intersection of the real part of 
the self-energy, and the short-dashed line is the linearized self¬ 
energy, which approximates the renormalized eigenvalue at 
the intersection of the x = y line. Lower: The corresponding 
spectral function. The position of the principal peak gives the 
quasiparticle energy. This narrow peak collects ~ 60% of the 
weight, and the rest of the charge forms a broad satellite peak 
below the bare eigenvalue. Note that the finite width of this 
peak is an artifact of the imaginary parameter used (0.1 eV). 


The quasiparticle broadening jx is defined as the 
half width of the spectral function at half of its max¬ 
imum, which means for a symmetrical quasiparticle 
peak that Aa(£a ± 7a) = ^a(£a)/2. Neglecting the fre¬ 
quency dependence of the self-energy near the quasipar¬ 
ticle energies, the broadening can be approximated as 
|7mE^^(£A)|, the imaginary part of the self-energy, which 
writes 


|7mEr(a;)| 


X! ^ 


X [{nv + fx')S{uj - Ey + uj^) 
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( 11 ) 


One recovers the static AHG expression for the ZPR 
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and the broadening by neglecting the phonon frequen¬ 
cies in the self-energy. Such approximation is made on 
the basis that the phonon frequencies and the quasiparti¬ 
cle corrections are much smaller than the typical energy 
differences with transition states that contribute to the 
self-energy. Under this assumption, the static Fan term 
reads 


static 
V Fan 


E 


1 I (V’aI IV^A') P 

y ^x - ^X' + sgn(eO) 


[2n,(T) + l]. 

( 12 ) 


Evaluating Eq. (12) instead of Eq. (4) requires much less 
computational efforts, especially if the Sternheimer equa¬ 
tion is used to eliminate the summation over high-energy 
electron bands. 

In this work, we adopt a semi-static approximation to 
compute the frequency-dependent self-energy. The terms 
of Eq. (4) are being computed explicitly, up to a certain 
band index and the contribution of the remaining 

bands above is treated statically with Eq. (12) using 
the Sternheimer equation method^^. For our materials, 
the bands above 5™^’' lie more than 20 eV above the 
states being corrected. Hence, the relative error on the 
self-energy due to the static treatment of these bands’ 
contributions is a few percent at most. 


Results and discussion 

We compute the quasiparticle structure of four crys¬ 
talline materials: diamond (C), boron nitride (BN) in 
the zinc-blende structure, magnesium oxide (MgO), and 
lithium fluoride (LiF) in the rock-salt structure. For all 
materials, we use a 8 x 8 x 8 k-point grid for the electronic 
wavefunctions and density, and a 32 x 32 x 32 q-point grid 
for the phonon modes sampling. 

The spectral functions at zero temperature are shown 
in Fig. 2 for the full band structure. A distinctive quasi¬ 
particle peak appears at the band edges, shifted from the 
bare eigenvalues, while in the regions of flat bands the 
spectral function is being diffused. The first conduction 
band of the indirect band gap materials (diamond, BN) 
exhibits a strong renormalization and a large broadening 
at F. This is due to the presence of other states in the 
Brillouin zone with close energies that are available for 
scattering. 

Another striking feature is the last valence band of LiF 
and MgO being completely diffused due to strong intra¬ 
band coupling. They show a narrow quasiparticle peak 
above the bare eigenvalue, and a broad satellite peak 
below. These features originate from the frequency de¬ 
pendence of both the real and the imaginary parts of the 
self-energy, as shown in Fig. 1 for the top of the valence 
bands of LiF. The satellite peaks could be observed in 
ARPES measurements, such as those performed on MgO 
by Tjeng et . However, a direct comparison with our 
results would require the full experimental spectra along 
the high-symmetry lines. 


Table I presents the real part of the self-energy for the 
states forming the optical band gap, namely the top of 
the valence bands (VB), and the first conduction band 
(CB) at F. At the bare eigenvalues, the self-energy shows 
little difference between the static and dynamical DFPT 
schemes, indicating that the phonon frequencies could 
be safely ignored in its real part. However, the frequency 
dependence produces an important renormalization fac¬ 
tor Z, ranging from 0.60 to 0.93 for the valence bands, 
and from 0.75 to 1.0 for the conduction bands. Thus, 
the dynamical effects tend to reduce the zero-point cor¬ 
rection, with respect to the static scheme. Comparing 
the linearized self-energy with the quasiparticle correc¬ 
tion obtained by solving Eq. (8) numerically or from the 
position of the principal peak of the spectral function, the 
linearization scheme proves to be a good approximation 
to both. 

The renormalization factor being larger than 1 indi¬ 
cates a breakdown of the quasiparticle picture. If the 
imaginary part of the self-energy is small, there is a well 
defined quasiparticle peak, and Z can be interpreted as 
the weight of that peak, which has to be smaller than 1. 
Otherwise, the spectral function is diffused, there is no 
such interpretation for Z and its value is unconstrained. 

Table H presents the quasiparticle broadening of 
the indirect-band gap materials computed with various 
schemes. The difference in the broadenings obtained 
from the static and the dynamical DFPT schemes is best 
understood with Eq. (11). Only the electronic states in 
a narrow energy range are available for scattering. The 
imaginary part of the self-energy is thus sensitive to the 
inclusion of phonon frequencies, since they affect the po¬ 
sitioning of this energy range. For the same reason, the 
broadening varies rapidly with frequency, which results 
in an important difference between the imaginary part at 
the bare eigenvalue and that at the renormalized eigen¬ 
value. Comparing these values with the actual width of 
the quasiparticle peak, we conclude that only the imagi¬ 
nary part of the self-energy evaluated at the renormalized 
energy is an accurate estimation of the broadening. 


II. ANHARMONIC EFFECTS 

The frozen-phonon method allows for a direct compu¬ 
tation of the electron-phonon self-energy within the adi¬ 
abatic approximation. We present here an extension of 
this method, which allows to explore anharmonic effects 
beyond the second-order perturbation theory of Allen, 
Heine and Cardona. 

Recalling the theory of the harmonic crystal, we write 
the total energy in a frozen-phonon configuration as 

E[^]=Eo + J2^zI, (13) 

where Eq is the equilibrium fixed-ions energy, z^, is a par¬ 
ticular phonon coordinate and z denotes the ensemble of 
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TABLE I. Zero-point renormalization (in eV) evaluated from the real part of the self-energy using 
a dynamical expression {dyn), at the bare eigenvalue (£°), at the renormalized eigenvalue (e), or 
the main quasiparticle peak from the bare eigenvalue (AA(£)). The unitless renormalization factor 
self-energy near the bare eigenvalue. See Apprendix A for the values of the rj parameter used. 


a static expression (stat), 
from the displacement of 
Z is used to linearize the 


\e°) 








AA(£) 


c 


BN 


MgO 


LiF 


VB 

CB 

Gap 

VB 

CB 

Gap 

VB 

CB 

Gap 

VB 

CB 

Gap 


0.134 

-0.238 

-0.372 

0.198 

-0.190 

-0.388 

0.197 

-0.153 

-0.350 

0.398 

-0.279 

-0.677 


0.126 

-0.240 

-0.366 

0.173 

-0.196 

-0.370 

0.198 

-0.143 

-0.341 

0.446 

-0.273 

-0.718 


0.931 

1.007 

0.823 

1.020 

0.734 

0.870 

0.596 

0.746 


0.118 

-0.242 

-0.359 

0.143 

- 0.200 

-0.343 

0.145 

-0.125 

-0.270 

0.266 

-0.204 

-0.469 


0.118 

-0.240 

-0.358 

0.147 

-0.197 

-0.344 

0.145 

-0.127 

-0.272 

0.254 

- 0.211 

-0.464 


0.118 

-0.247 

-0.365 

0.147 

-0.208 

-0.355 

0.147 

-0.127 

-0.274 

0.256 

- 0.211 

-0.467 



BN band structure (eV) 



LiF band structure (eV) 



MgO band structure (eV) 



FIG. 2. Spectral functions summed over the bands at each k-point of the Brillouin zone (arbitrary units). The green lines are 
the DFT band structure, in eV. When a quasiparticle peak is visible in the spectral function, the renormalization is infered 
from the difference of the position of that peak with the bare band structure. In the regions of flat bands, the band structure 
is being completely diffused. A satellite peak is seen below the last valence band of LiF and MgO. 
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TABLE II. Quasiparticle broadening (in eV) evaluated from 
the imaginary part of the self-energy using a static expression 
(stat), a dynamical expression {dyn), at the bare eigenvalue 
(e®), at the renormalized eigenvalue (e), or from the width of 
the main quasiparticle peak at half of its maximum ( 7 ). See 
Apprendix A for the values of the rj parameter used. 



pmE'’*“*(£“)| 

pmE‘^J'"(£°)| 

pmE‘^^’"(e)| 

7 

C CB 

0.178 

0.164 

0.140 

0.138 

BN CB 

0.246 

0.226 

0.200 

0.196 


FIG. 3. Some of the high-order electron-phonon coupling dia¬ 
grams which contribute to the self-energy. Each vertex with n 
phonon branches is associated with a n*^-order derivative of 
the one-particle Hamiltonian. The independent-phonon ap¬ 
proximation presented in the text retains only the diagrams 
formed with multiple Interactions involving the same phonon 
mode. 

all of these coordinates. Taking the lattice dynamics into 
account, the phonon eigenstates are those of the decou¬ 
pled harmonic oscillators: 

V 

where n denotes the ensemble of all the phonon occupa¬ 
tion numbers. The total energy in this state is 

E[n] = Eq+ ^uj„[n^ + ^]. (15) 

The expression for a particular eigenvalue at finite tem¬ 
perature is given by the derivative oi E = E — TS, the 
Helmholtz free energy with respect to an electronic oc¬ 
cupation number f\, which reduces to"*’^ 

= = + (16) 

This expression should be compared with the electron- 
phonon self-energy in the adiabatic approximation 
(Eq. (5) and (12).) The individual phonon contributions 
to the self-energy are proportional to du/^ldfx, which 
we call the electron-phonon coupling energies (EPCE). 
Using Brook’s theorem^^, the EPCEs are given by the 
second-order derivatives of an eigenvalue with respect to 
a phonon coordinate: 

duji, d£\ 1 r 

dfx dn^ 2io„ dzl 


where ex [z^,] is an electronic energy computed with all 
atoms displaced by a length z^, along the phonon dis¬ 
placement vector [/T. This expression does not rely 
on the rigid-ion approximation, but requires a super¬ 
cell calculation to account for the phonon wavevector. 
Within the validity of the rigid-ion approximation, it 
should reproduce the results of the static DEPT scheme. 
Both of these frameworks are developed within the har¬ 
monic approximation, since the total energy and the elec¬ 
tronic eigenvalues are expanded up to second order in the 
phonon perturbations. Equivalently, the harmonic ap¬ 
proximation can be defined as the assumption that the 
electronic eigenvalues vary quadratically with a phonon 
coordinate z^. 

In order to relax the harmonic approximation on the 
electronic energies, we cast the free energy E = ksT In Z 
in terms of the canonical partition function Z = Tr , 
which is a trace over both the electronic and the atomic 
degrees of freedom. Resolving the trace over the electron 
coordinates, the expression for a temperature-dependent 
eigenvalue reads 

r p-PE[z,] 

£\{T) = J dz—^—ex[z], (18) 

where Zi = J dz is the partition function of 

the atoms only, and £x[z] is the eigenvalue computed 
in some frozen-phonon configuration. This formulation 
is reminiscent of the path-integral molecular dynamics 
approach^®’^®, with the difference that a configuration 
is specified in terms of phonon coordinates rather than 
atomic positions in real space. It retains the adiabatic 
approximation, since the atomic motion does not induce 
electronic transitions, but leaves the electrons in their 
evolving states. 

We now use the crystal phonon structure and perform 
the harmonic approximation on the total energy only, 
writing 

g-/3£[n] r 

£\{T) = Y^ J dz\x’^{z)\ ex[z]. (19) 

We are assuming here that the set of phonon wavefunc- 
tions (zv) and frequencies lov computed from second- 
order perturbation theory are good eigenfunctions of the 
system. That is to say that the total energy is quadratic 
along the computed phonon modes even if the eigenval¬ 
ues are not. Equation (19) now includes all high-order 
diagrams that may contribute to the self-energy, such as 
those depicted in Fig. 3. 

Finally, we perform the independent phonon approxi¬ 
mation and write 

£A(T) = eUEE"”‘' / dz.\x--'iz,)\\ex[z.]-el), 

u riu 

( 20 ) 

where s”'" = In doing so we 

disregard the cross-terms contributions between different 
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phonons modes. This ansatz restricts the high-order di¬ 
agrams to those containing a single phonon mode, which 
may interact multiple times with the electrons. These 
additional diagrams come from the anharmonicity of the 
eigenvalues appearing in the integrant of Eq. (20), as il¬ 
lustrated in Fig. 4. One can verify that if the eigenvalues 
vary quadratically with the phonon displacements, then 
Eq. (16) is recovered. Otherwise, Eq. (20) defines ef¬ 
fective EPCEs for each phonon mode which include the 
anharmonic effects. 


Results and discussion 

We compute the EPCEs by frozen-phonon calcula¬ 
tions, using the phonon displacement vectors obtained 
from DEPT. For the harmonic approximation, Eq. (17) 
is evaluated with atomic displacements of about 10“^A, 
while the anharmonic effects are included by evaluating 
Eq. (20) with 20 displacements up to ~ O.SA, which cor¬ 
responds to about 4 units of a typical phonon average 
displacement l/y^av. 

The EPCEs are shown in Fig. 5 through the Brillouin 
zone of diamond. The spiky structure of the EPCEs of 
the first conduction band at P results from this state 
not being at the bottom of the band. Consequently, 
when a phonon wavevector connects the state at P to an 
other state with close energy, a divergence occurs in the 
EPCEs. A divergence also occurs for phonon wavevec- 
tors near P, for both the VB and CB states, but these 
divergences integrate to a finite value when the density of 
phonon modes is taken into account. The EPCEs com¬ 
puted with the frozen-phonon method in the harmonic 
approximation are in close agreement with the DEPT re¬ 
sults, indicating that the rigid-ion approximation holds. 
However, when the full dependence of the eigenvalues on 
the phonon displacements is taken into account, the an¬ 
harmonicity of the eigenvalues tend to reduce the EPCEs, 
with respect to the harmonic approximation. 

This is exemplified on Fig. 4 with the mode 174 (the 
fourth mode with wavevectors = {L + X)/2). In the 
harmonic approximation, it contributes —869 meV to the 
CB EPCE at this q-point. The eigenvalue however de¬ 
parts from the quadratic behavior with the phonon dis¬ 
placement, reducing the coupling energy to —383 meV. 
On the other hand, the total energy follows closely the 
quadratic curve, indicating that this displacement is a 
genuine phonon mode. This tendency is observed near all 
divergent points of the Brillouin zone and near the zone 
center. The second-order perturbation theory is thus in¬ 
sufficient to treat the effect of those strongly coupling 
modes on the electronic states. 

Table III reports the ZPR computed on a 4 x 4 x 4 q- 
point grid with the various static schemes. Again, the 
total ZPR obtained with the harmonic frozen-phonon 
method and with DEPT are in good agreement. The 
discrepancies can be attributed to the rigid-ion approxi¬ 
mation. When anharmonic effects are included, the total 


TABLE III. Zero-point renormalization of the band gap (in 
eV) within the adiabatic approximation, obtained whith the 
static DFPT, with the frozen-phonon method in the harmonic 
approximation (FPH), and with the frozen-phonon method 
including anharmonic effects (FPA). 



Static DFPT 

FPH 

FPA 

c 

VB 

0.115 

0.119 

0.107 


CB 

-0.320 

-0.321 

-0.214 


Gap 

-0.436 

-0.439 

-0.320 

BN 

VB 

0.120 

0.133 

0.108 


CB 

-0.193 

-0.198 

-0.154 


Gap 

-0.313 

-0.331 

-0.262 

MgO 

VB 

0.110 

0.118 

0.070 


CB 

-0.081 

-0.078 

-0.084 


Gap 

-0.191 

-0.196 

-0.154 

LiF 

VB 

0.445 

0.431 

0.168 


CB 

-0.130 

-0.122 

-0.113 


Gap 

-0.575 

-0.553 

-0.281 


renormalization of the electronic energies is typically re¬ 
duced compared to the harmonic approximation. For the 
indirect band gap materials, the renormalization of the 
CB state is largely affected by the anharmonic effects, 
since they receive an important contribution from those 
strongly coupling modes at the Brillouin zone boundaries 
which are being attenuated. The states at the band edges 
are being affected to various extends. The valence band 
of LiF, which is especially flat, shows a strong anhar¬ 
monicity in the ZPR coming from the modes near P, 
reducing the ZPR by about 60%. In contrast, the con¬ 
duction band of MgO, which is very dispersive, is only 
slightly affected by these effects. 

Our results are obtained on a coarse q-point grid, 
limited by the scaling of the frozen-phonon method. 
Whether the same conclusions apply to a converged q- 
point grid depends on the relative importance of strongly 
coupling modes, since they are responsible for anhar¬ 
monic effects. For the states lying at the top or at 
the bottom of their respective band, the ZPR increases 
monotonically with the number of q-points on a regu¬ 
lar grid. This is because the region of strongly coupling 
phonon modes near T gains in importance. Thus, the an¬ 
harmonic effects for these state are expected to grow as 
a finer q-point sampling is achieved. For the states that 
are not at the extrema of their band, the convergence of 
the ZPR with q-points sampling is non-monotonic. In 
these cases, we cannot make quantitative predictions for 
the anharmonic effects on a converged q-point sampling. 
Nevertheless, the presence of strongly coupling modes 
near the Brillouin zone boundaries suggests an impor¬ 
tant anharmonic contribution to the ZPR. 
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FIG. 4. Dependence of the CB eigenvalue (green, circles) and 
the total energy (blue, squares) on a phonon displacement for 
the mode D 4 of diamond. The circles and squares are the 
actual frozen-phonon calculations, the solid lines correspond 
to the harmonic approximation, and the filled curve is the 
phonon wavefunction. 


III. CONCLUSION 

The dynamical DFPT scheme allowed us to compute 
the frequency-dependent electron-phonon coupling self¬ 
energy. Our calculations yield a renormalization factor 
ranging from 0.6 to 1.0. This renormalization factor is 
important to obtain correct quasiparticle energies, but 
has been laregly overlooked in the literature. 

The spectral function reveals distinctive features of the 
quasiparticle band structure. In the indirect band gap 
materials (diamond and BN), the conduction bands un¬ 
dergo intra-band scattering processes, which broaden the 
spectral function at T. In the direct band gap materials 
with flat valence bands (LiF and MgO), these processes 
even generate satellite peaks below the valence bands. 

The broadening can be obtained from the imaginary 
part of the self-energy, but one has to use a dynamical 
theory to do so. Not only are the phonon frequencies 
necessary to impose energy conservation in the scattering 
process, but the imaginary part of the self-energy must 
be evaluated at the renormalized eigenvalues, in order to 
compute properly the quasiparticle broadening. 

Finally, we explored anharmonic effects using frozen- 
phonon calculations. The anharmonicity in the eigen¬ 
value dependence on the atomic displacements occurs 
even if the phonon modes are correctly described by the 
second-order perturbation theory. This effect tend to de¬ 
crease the contribution of the strongly coupling phonon 
modes, reducing the ZPR of certain states by as much as 
60% with respect to the static AHC theory. 

Our results indicate that high-order electron-phonon 
coupling terms bring an important contribution to the 
self-energy and the ZPR. Our methodology however in¬ 
cludes a partial summation of the high-order terms and 
treats the perturbations statically. A theory that would 
include all high-order terms in a dynamical way cannot 
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FIG. 5. Upper: The band structure of diamond, in eV. The 
dashed lines show the energies of the VB and CB states. Mid¬ 
dle and lower: Electron-phonon coupling energies (EPCE), 
in eV for the CB state (middle) and the VB state (lower), 
computed with various methods. The blue line is the DFPT 
calculation, the yellow discs are the frozen-phonon method 
in the harmonic approximation, and the green triangles are 
the frozen-phonon method including anharmonic effects. A 
divergence is observed in the EPCEs of the CB state when 
a phonon wavevector couples this state to an other one with 
close energy, while the EPCE of both VB and CB states show 
a broad diverging peak at the center of the Brillouin zone. 


be tested at present, but could be eventually addressed 
with quantum Monte-Carlo approaches. 
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Appendix A: Imaginary parameter and convergence 
properties 

In order to compute the ZPR, one has to sample the 
phonon wavevectors in the Brillouin zone, either through 
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FIG. 6. Upper: Convergence of the self-energy at the renor¬ 
malized eigenvalue for the CB state of diamond, as a func¬ 
tion of q-points spacing, with various imaginary parameters. 
Lower: Frequency dependence of the self-energy for the CB 
state of diamond, with various imaginary parameters. The 
solid lines are obtained on a 32 x 32 x 32 q-point grid, and 
the dashed lines are obtained on a 24 x 24 x 24 q-point grid, 
corresponding to the two left-most points on the upper pannel 
graph. 


a regular mesh, or with a random set of q-points. At the 
same time, one has to select a value for the parameter rj 
giving an imaginary part to the self-energy. The choice of 


this parameter should in facts be addressed in conjunc¬ 
tion with the q-points sampling. When the static DFPT 
method is used, the numerical value assigned to r] is usu¬ 
ally on the order of typical phonon frequencies 0.1 eV) 
to account for their omission. Otherwise, if a too small 
value of r] is used, it is not even clear that the self-energy 
will converge to a finite value^^. In a dynamical scheme, 
one should in principle aim for a vanishing value for rj. 
While it is expected that the self-energy elements should 
converge to a finite value as 77 —>■ 0, tuning the value of 
77 conveniently eases the convergence with the number of 
q-points, as shown in Fig. 6 . 

Moreover, using a too small value of 77 can compromise 
the frequency dependence of the self-energy, as shown on 
Fig. 6 for the CB state of diamond. Even for the most 
converged q-point grid, the self-energy computed with 
77 = 0.1 eV shows rapid variations with uj. These vari¬ 
ations are even larger for a smaller q-point grid, and in 
these cases, the solution of w = E(w -|-e°) could certainly 
not be estimated by linearizing the self-energy near the 
bare eigenenergy. The self-energy becomes a perfectly 
smooth function of uj when 77 = 0.4 eV. 

We use the following criterion to determine the value 
of 77 . Consider the contribution of two neighboring q- 
points q and q' to the self-energy of the electronic state 
hn. The contribution of a particular electron band n 
and phonon branch m will have terms proportional to 
[(w - £k+qn ± + {uj - Ek+q'n ± ^q'm + , 

assuming that the matrix elements in the numerator of 
the self-energy does not change between q and q'. If 
the value of 77 is vanishingly small, the spectral func¬ 
tion will exhibit distinct peaks at a; = Ek+qn ‘^qm 
and Lu = Ek-fq'n ^ <^q'm, which will be an artifact of 
the q-points sampling. The separation of those peaks 
comes mainly from the dispersion of the electronic ener¬ 
gies, which is more important than that of the phonon 
frequencies. Simple analysis shows that these peaks can 
be made undistinguishable by setting 77 = -x/SAe /2 where 
Ae = Ek+qn ~ ^k-i-q'n' H^nce, for a given q-point mesh, 
we compute the largest Ae between neighboring q-points, 
within the bands being corrected, and use it to deduce 77 . 
The values of 77 obtained for our most converged q-point 
grid are: 0.2 eV for the VB state of diamond, 0.4 eV for 
the CB states of diamond and BN, and 0.1 eV for all the 
other VB and CB states. We verified that the broaden¬ 
ing of the CB states of diamond and BN was insensitive 
to the choice of this parameter. 
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